!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine source_term (mx_new, my_new, mz_new, bx_new, by_new, bz_new, en_new, az_new, po_new,                                    &
      ro_old, mx_old, my_old, bx_old, by_old, bz_old, po_old, gx, gy,                                                              &
      bx_resis, by_resis, bz_resis, en_resis, az_resis,                                                                            &
      ch2, cp2, dx, dy, dt, nx, ny, div_clean_flag)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny, div_clean_flag

   real(8), intent(in) :: dx, dy, dt, ch2, cp2

   real(8), dimension(nx, ny), intent(inout) :: mx_new, my_new, mz_new, bx_new, by_new, bz_new, en_new, az_new, po_new
   real(8), dimension(nx, ny), intent(in) :: ro_old, mx_old, my_old, bx_old, by_old, bz_old, po_old, gx, gy
   real(8), dimension(nx, ny), intent(in) :: bx_resis, by_resis, bz_resis, en_resis, az_resis !, te_cond_x, te_cond_y

   integer(4) :: i, j, ip1, im1, jp1, jm1
   real(8) :: div_b, s_mx, s_my, s_mz, s_bx, s_by, s_bz, s_en, s_az, dx2, dy2

!----------------------------------------------------------------------------------------------------------------------------------|
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0

   do j = 2, ny - 1
      jp1 = j + 1
      jm1 = j - 1
      do i = 2, nx - 1
         ip1 = i + 1
         im1 = i - 1

         div_b = ((bx_old(ip1, j) - bx_old(im1, j)) / dx2 + (by_old(i, jp1) - by_old(i, jm1)) / dy2) * div_clean_flag
         s_mx = -div_b * bx_old(i, j) + ro_old(i, j) * gx(i, j)
         s_my = -div_b * by_old(i, j) + ro_old(i, j) * gy(i, j)
         s_mz = -div_b * bz_old(i, j)
         s_bx = bx_resis(i, j)
         s_by = by_resis(i, j)
         s_bz = bz_resis(i, j)
         s_en = -(po_old(ip1, j) - po_old(im1, j)) / dx2 * bx_old(i, j) * div_clean_flag -                                         &
                 (po_old(i, jp1) - po_old(i, jm1)) / dy2 * by_old(i, j) * div_clean_flag +                                         &
!                 (te_cond_x(ip1, j) - te_cond_x(im1, j)) / dx2 + (te_cond_y(i, jp1) - te_cond_y(i, jm1)) / dy2 +                   &
                 mx_old(i, j) * gx(i, j) + my_old(i, j) * gy(i, j) + en_resis(i, j)
         s_az = (mx_old(i, j) * by_old(i, j) - my_old(i, j) * bx_old(i, j)) / ro_old(i, j) + az_resis(i, j)

         mx_new(i, j) = mx_new(i, j) + s_mx * dt
         my_new(i, j) = my_new(i, j) + s_my * dt
         mz_new(i, j) = mz_new(i, j) + s_mz * dt
         bx_new(i, j) = bx_new(i, j) + s_bx * dt
         by_new(i, j) = by_new(i, j) + s_by * dt
         bz_new(i, j) = bz_new(i, j) + s_bz * dt
         en_new(i, j) = en_new(i, j) + s_en * dt
         po_new(i, j) = po_new(i, j) * exp (-dt * ch2 / cp2 * div_clean_flag)
         az_new(i, j) = az_new(i, j) + s_az * dt
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine source_term
